← LAB ◊ ↑ TOC ↑

Name, Vorname: Taubenberger, Samuel

LAB 2: Filterung und DFT als Näherung für spektrale Analyse

Inhalt

0. Allgemeines

1. FIR und IIR Filter

2. Schätzung der Fourierreihe durch FFT

3. Schätzung des Fourierspektrums von Sinustönen mit überlagertem Rauschen durch FFT

4. Fourierspektrum eines Rechteckpulses

0. Allgemeine Hinweise

• <SHIFT>-<RETURN> führt eine Codezelle aus und rendert eine Textzelle.
• In Markdown sind Leerzeilen wichtig zum Trennen von Abschnitten!
• Sie können LaTeX-Code zwischen \$ ... \$ einschließen.
• Kontexthilfe zu Funktionen etc. bekommen Sie über <SHIFT>-<TAB>

Nach dem Praktikumsversuch exportieren Sie das Notebook mit Textantworten, Codezellen und Plots als HTML (File -> Export Notebook As ... -> Export Notebook to HTML) und reichen es in Moodle ein.

---

import os, sys
module_path = os.path.abspath(os.path.join('..')) # append directory one level up to import path
if module_path not in sys.path:  # ... if it hasn't been appended already
    sys.path.append(module_path)
import dsp_fpga_lib as dsp
dsp.versions() # print versions 
    
%matplotlib inline
import matplotlib.pyplot as plt
size = {"figsize":(12,5)} # Plotgröße in Inch

import numpy as np
import scipy.signal as sig
import wave

from IPython.display import Audio, display
#-----------------------------------------------------------------------------
def wav2np(filename):
    """ Read the wav-file and convert it to a one or two-dimensional numpy array,
        depending on the number of channels. 
        Properties of the WAV-file are stored as function attributes (evil)
    """
    wf = wave.open(filename,'rb')
    wav2np.N_CH = wf.getnchannels() # number of channels
    wav2np.W    = wf.getsampwidth() # wordlength per sample in bytes
    wav2np.N_FR = wf.getnframes()   # number of frames
    wav2np.f_S  = wf.getframerate() # sample (frame) rate 
    
    print("{0} channels with {1} frames of {2} bytes and f_S = {3} Hz.".format(wav2np.N_CH, wav2np.N_FR, wav2np.W, wav2np.f_S))
    
    if wav2np.W == 2:
        samples_in = np.frombuffer(wf.readframes(-1), dtype=np.int16) # read wav data as 16 bit integers, R and L samples are interleaved
    elif wav2np.W == 1:
        samples_in = np.frombuffer(wf.readframes(-1), dtype=np.int8) # read wav data as 8 bit integers, R and L samples are interleaved
    else:
         raise TypeError("Unknown data format: {0} bytes".format(wav2np.W))   

    samples = np.array([samples_in[idx::wav2np.N_CH] for idx in range(wav2np.N_CH)], dtype=np.int32) # deinterleave channels to numpy array N_CHAN x N_FRAMES
    return samples

Python version: 3.12.0
Numpy: 1.26.0
Scipy: 1.11.3
Matplotlib: 3.8.0 module://matplotlib_inline.backend_inline

1. FIR und IIR Filter

In diesem Abschnitt analysieren Sie die beiden Filter aus der folgenden Abbildung:

Filter 1		Filter 2

VORBEREITUNG:

• Woran erkennen Sie, welches davon ein IIR- und welches ein FIR Filter ist?

• Berechnen Sie die Differenzengleichung der beiden Filter.

• Berechnen Sie die Systemantwort der beiden Filter. Können Sie die Impulsantwort beider Filter einfach aufschreiben?

• Berechnen Sie Pole und Nullstellen beider Systemfunktionen und skizzieren Sie (auf einem Blatt Papier) den P/N - Plan.

• Schätzen Sie aus dem P/N Plan ab, um welchen Filtertyp es sich handelt. Bei welcher Frequenz liegt das Maximum des Betragsgangs? Welchen Wert hat es?

---

1.1 Simulation

Für die Simulation müssen wir zuerst die Filter als Koeffizientenvektoren a und b darstellen (siehe folgendes Bild und Abschnitt 3 von LAB 1).

Als "Codesteinbruch" können Sie wieder das Notebook 02_LTF/LTF-Filter_properties.ipynb verwenden.

Alternativ verwenden Sie pyfda: Im Tab "b,a" importieren Sie Koeffizienten aus einem CSV-File oder (nach Auswahl von "Clipboard" in den Einstellungen) direkt aus der Zwischenablage. Die Koeffizienten b des nicht-rekursiven Teils geben Sie hierfür getrennt durch Kommata an, optional in einer zweiten Zeile die Koeffizienten a des rekursiven Teils.

SIMULATION:

• Stellen Sie Filter 1 und Filter 2 (s.o.) als Koeffizientenvektoren a und b dar. Informationen dazu finden Sie im Abschnitt 3 von LAB 1 und in obiger Abbildung. Beachten Sie:

  • Die Reihenfolge ist $a = [1, a_1, a_2, ...]$ und $b = [b_0, b_1, \ldots]$. Warum ist immer $a_0 = 1$? Warum dreht sich das Vorzeichen der rekursiven Koeffizenten herum? Antwort: Kommt durch Ausklammern von Y(z). Da wir den rekursive Anteil zur Aufstellung der Übertragungsfunktion immer auf die andere Seite bringen.
  • "Fehlende" Koeffizienten müssen als 0 eingetragen werden.
  • Die imaginäre Einheit ist 1j, dementsprechend werden imaginäre Zahlen z.B. als 0.3j repräsentiert.

• Testen Sie mit dem P/N Diagramm, ob Rechnung und Simulation zusammen passen.

• Lassen Sie sich die Impulsantwort anzeigen, für das FIR-Filter ist das relativ einfach ;-), beim IIR-Filter benötigen Sie dsp.impz()

• Plotten Sie Betragsfrequenzgang, Phasengang und Gruppenlaufzeit.

1.1.1 P/N Diagramm

a_FIR = [1,0,0,0]
b_FIR = [0,0.5,1,0.5]
a_IIR = [1,0,0.9] #see function definition of e.g. dsp.impz()
b_IIR = [0,0,1]
#---------------
fig, (ax1,ax2) = plt.subplots(nrows=1, ncols=2,**size)
ax1.set_title('P/N Diagramm: FIR')
dsp.zplane(b_FIR,a_FIR, plt_ax=ax1);
print("Nullstellen FIR: {0}".format(np.roots(b_FIR)))

ax2.set_title('P/N Diagramm: IIR')
dsp.zplane(b_IIR,a_IIR, plt_ax=ax2);
print("Nullstellen IIR: {0}".format(np.roots(b_IIR)))
if type(a_IIR) in {list, np.ndarray} and len(a_IIR) > 1:
    print("Polstellen IIR:  {0}\n".format(np.roots(a_IIR)))

Nullstellen FIR: [-1. -1.]
Nullstellen IIR: []
Polstellen IIR:  [-0.+0.9486833j  0.-0.9486833j]

1.1.2 Impulsantwort

Der Befehl h,t = dsp.impz(b,a) berechnet die Impulsantwort, die Sie dann plotten können, am Besten als stem-Plot.

Warum ist beim IIR-Filter jeder zweite Impuls Null? Antwort: Weil beim Rückführungszweig die Komponente mit z^-1 fehlt und nur z^-2 vorhanden ist.

fig, (ax1,ax2) = plt.subplots(nrows=2, ncols=1,**size)
ax1.set_title('Impulsantwort FIR')
h_FIR,t_FIR = dsp.impz(b_FIR,a_FIR);
h_IIR,t_IIR = dsp.impz(b_IIR,a_IIR);
ax1.stem(t_FIR, h_FIR)
ax2.set_title('Impulsantwort IIR')
ax2.stem(t_IIR, h_IIR)
fig.set_tight_layout(True)

1.1.3 Betrags- und Phasengang sowie Gruppenlaufzeit

Aus dem komplexen Frequenzgang omega, H = sig.freqz(b,a) ermitteln Sie mit np.abs() und np.angle() Betrags- und Phasengang. Defaultmäßig werden 512 Frequenzpunkte zwischen 0 und $f_S/2$ bestimmt und als normierte Kreisfrequenz zurückgegeben.

Die Gruppenlaufzeit ermitteln Sie mit w, H = sig.group_delay((b,a), omega). Mit omega übergeben Sie die Frequenzen, an denen die Gruppenlaufzeit berechnet werden soll (z.B. die, die Sie aus der Berechnung des komplexen Frequenzgangs erhalten haben).

Der Plot der Gruppenlaufzeit sieht u.U. etwas seltsam aus, mit set_ylim([y_min, y_max]) passen Sie die Grenzen an.

fig, ax = plt.subplots(nrows=3, ncols=2, figsize=(12,8))
w_FIR, H_FIR = sig.freqz(b_FIR, a_FIR)
_, tau_FIR = sig.group_delay((b_FIR, a_FIR), w_FIR)
w_IIR, H_IIR = sig.freqz(b_IIR, a_IIR)
_, tau_IIR = sig.group_delay((b_IIR, a_IIR), w_IIR)
ax[0][0].set_title("FIR - Betragsgang")
ax[0][0].plot(w_FIR / (2*np.pi), np.abs(H_FIR))
ax[0][0].set_ylabel(r"$|H| \;\rightarrow$")

ax[0][1].set_title("FIR - Phasengang")
ax[0][1].plot(w_FIR / (2*np.pi), np.angle(H_FIR))
ax[0][1].set_ylabel(r"$phi \;\rightarrow$")

ax[1][0].set_title("FIR - Gruppenlaufzeit")
ax[1][0].plot(w_FIR/ (2*np.pi), tau_FIR)
ax[1][0].set_ylabel(r"$tau \;\rightarrow$")

ax[1][1].set_title("IIR - Betragsgang")
ax[1][1].plot(w_IIR / (2*np.pi), np.abs(H_IIR))
ax[1][1].set_ylabel(r"$|H| \;\rightarrow$")

ax[2][0].set_title("IIR - Phasengang")
ax[2][0].plot(w_IIR / (2*np.pi), np.angle(H_IIR))
ax[2][0].set_ylabel(r"$phi \;\rightarrow$")

ax[2][1].set_title("IIR - Gruppenlaufzeit")
ax[2][1].plot(w_IIR/ (2*np.pi), tau_IIR)
ax[2][1].set_ylabel(r"$tau \;\rightarrow$")


ax[2][1].set_xlabel(r"$F \;\rightarrow$")   
fig.set_tight_layout(True)

1.2 Anhören

Im folgenden filtern wir Sprach- oder Rauschsignale mit unseren beiden Filtern und hören und schauen uns das Resultat an. Eine FIR-Filterung könnten Sie wieder mit convolve(x,h) durchführen (nur für eindimensionale Arrays). Warum funktioniert bei IIR-Filtern die convolve-Methode nicht? Antwort: Würde eine unendlich ausgedehnte Impulsantwort hervorrufen.

Für IIR und FIR-Filter können Sie die Routine y = sig.lfilter(b,a,x) verwenden (funktioniert auch mit zweidimensionalen Arrays).

SIMULATION:

• Filtern Sie ein Rauschsignal oder einen WAV-File aus dem Unterordner medien mit dem FIR- und dem IIR-Filter. Rauschen erzeugen Sie wieder z.B. mit x_n = np.random.randn(16000), WAV-Files wandeln Sie mit der Hilfsfunktion wav2array(filename) in ein ein- oder zweidimensionales Array um.

• Hören Sie sich die Wirkung der beiden Filter an mit der Audio Klasse aus dem IPython.display - Modul:

  display(Audio((data=None, rate=None), data kann dabei ein ein- oder zweidimensionales numpy-Array oder Liste sein, ein Filename oder auch eine URL. Der Parameter rate definiert die Abtastrate.

• Schauen Sie sich einen Ausschnitt des ursprünglichen und des gefilterten Signals im gleichen Plot an und vergleichen Sie.

"""Lab1: 
x_n = np.squeeze(np.random.randn(50000,1))
h = np.ones(16) # MA-filter mit Länge L - setzen Sie den passenden Zahlenwert ein
y = np.convolve(x_n, h)
len(y) #500+16-1
display(Audio(data=y.T, rate=16000))"""

x_n = np.random.randn(16000)
x_w = wav2np("../medien/87778__marcgascon7__vocals.wav")



print("Rauschen original:")
display(Audio(data=x_n, rate=wav2np.f_S))

print("WAV (original):")
display(Audio(data=x_w, rate=wav2np.f_S))

print("Gefilteretes rauschen:")
print(np.shape(x_w))
y_FIR_R = sig.lfilter(b_FIR, a_FIR, x_n)
y_IIR_R = sig.lfilter(b_IIR, a_IIR, x_n)

print("FIR:")
display(Audio(data=y_FIR_R, rate=wav2np.f_S))
print("IIR:")
display(Audio(data=y_IIR_R, rate=wav2np.f_S))

print("Gefiltertes WAV:")
print(np.shape(x_w))
y_FIR_W = sig.lfilter(b_FIR, a_FIR, x_w)
y_IIR_W = sig.lfilter(b_IIR, a_IIR, x_w)

print("FIR:")
display(Audio(data=y_FIR_W, rate=wav2np.f_S))
print("IIR:")
display(Audio(data=y_IIR_W, rate=wav2np.f_S))

2 channels with 349952 frames of 2 bytes and f_S = 44100 Hz.
Rauschen original:

WAV (original):

Gefilteretes rauschen:
(2, 349952)
FIR:

IIR:

Gefiltertes WAV:
(2, 349952)
FIR:

IIR:

fig, ax = plt.subplots(**size)
l = 4000
start = 20000
print(np.arange(l))
print(len(y_IIR_W))
ax.plot(np.arange(l), x_w[0][start:start+l], label="x_w")
ax.plot(np.arange(l), y_IIR_W[0][start:start+l], label="y_IIR")
ax.legend()

[   0    1    2 ... 3997 3998 3999]
16000

<matplotlib.legend.Legend at 0x2046b6905f0>

2. Fourierreihe einer Rechteckpulsfolge

In diesem Versuchsteil bestimmen wir die Fourierkoeffizienten einer rechteckförmigen zeitkontinuierlichen Pulsfolge mit $T_1 = 1$ ms. Der Duty Cycle $\alpha = 0.25$ und Amplitude $A=1000$ mV, die Pulsbreite ist also $\Delta T = T_1/4$ (siehe nächster Plot).

fig, ax = plt.subplots(**size)
t = np.arange(-0.2e-3, 4e-3, 1/640e3) # pseudo-analoger Zeitvektor in ms (mit 640 kHz so fein abgetastet, dass der Unterschied nicht auffällt)
y = 500*sig.square(t*1e3*2*np.pi+np.pi/4, duty=0.25) + 500
ax.plot(t, y, 'r', lw=2)
ax.axvline(0, ls='-',color='k')
ax.set_xlabel(r"$t \; / \;\mathrm{s} \;\rightarrow$")
ax.set_ylabel(r"$y(t)\; / \;\mathrm{mV} \;\rightarrow$");

2.1 Berechnung der Koeffizienten

Ein mit $T_1$ periodisches Signal $y(t)$ lässt sich mit Hilfe der Fourierreihe für reellwertige Signale (Index $k \in \mathbb{N}_0$) darstellen:

$$ y(t)= {a_0} + 2\sum_{k=1}^\infty a_k \cos 2 \pi k f_1 t + b_k \sin 2 \pi k f_1 t \; k \in \mathbb{N} $$

Zunächst berechnen wir die Koeffizienten $a_k$, $b_k$ für $k = 0, 1, \ldots, 11$ mit Hilfe der Formel für die Fourierreihe von reellwertigen Zeitsignalen. Da das Signal achsensymmetrisch ist, sind die imaginärwertigen Koeffizienten $b_k=0$.

$$ \begin{align}{a}_k &= \frac 1 {T_1} \int_{-T_1/2}^{T_1/2}y(t)\cos(2 k \pi f_1 t) \, \text{d}t = \frac {A} {T_1} \int_{-\alpha T_1/2}^{\alpha T_1/2}\cos(2 k \pi f_1 t) \, \text{d}t \\ &= \left. \frac {A} {T_1 \cdot 2 k \pi f_1} \sin(2 k \pi f_1 t)\, \right|_{-\alpha T_1/2}^{\alpha T_1/2} = {\frac {A} { k \pi} \sin ( k \alpha \pi) }= { {A \alpha}\, \mathrm{si} ( k \alpha \pi) } \end{align}$$

SIMULATION:

Schreiben Sie ein Skript, das die ersten 11 Koeffizienten der Rechteckfunktion exakt berechnet. Welchen Vorteil hat es, die sinc(x) Funktion zu verwenden anstelle von sin(x)/x? Aber Achtung: In Numpy (und in Matlab) ist die sinc-Funktion definiert als $\mathrm{sinc}(x) = sin(\pi x)/\pi x$, das $\pi$ im Argument müssen Sie daher weglassen. Drucken Sie die Ergebnisse übersichtlich in eine Tabelle aus (s.u.)

2.1.1 Formatierte Ausdrucke in Python

Einen formatierten Ausdruck erhalten Sie z.B. mit print("\n{0:7.2f} | ".format(Y[i]), end="") (der Teil in der geschweiften Klammer wird ersetzt durch Y[i] und formatiert mit insgesamt 7 Stellen, 2 Nachkommastellen, keinen Zeilenumbruch). Eingebettet in eine for Schleife bekommen Sie so schnell eine übersichtliche Tabelle.

• "\n" erzeugt einen Zeilenumbruch in einem String

• print(a, end="") ersetzt den Zeilenumbruch am Ende des Printbefehls durch ein anderes Zeichen (oder nichts wie hier)

• Die str.format() Anweisung ersetzt Ausdrücke in geschweiften Klammern im String durch die Argumente in der .format() Anweisung, versuchen Sie es selbst mit "{0} ich heiße und bin {1} Jahre alt".format("Yoda", 877).

• Der Ausdruck in der geschweiften Klammer kann durch Formatierungsanweisung ergänzt werden, {0:7.2f} formatiert das Argument 0 als Float mit 2 Nachkommastellen und insgesamt mindestens 7 Stellen. {0:>7} lässt ebenfalls Platz für mindestens 7 Stellen oder Zeichen und richtet den Ausdruck rechtsbündig aus.

  Auf https://www.python-kurs.eu/python3_formatierte_ausgabe.php finden Sie eine übersichtliche Grafik hierzu.

k=np.arange(11)
alpha = 4/16; A = 1000
Y0 =  A*alpha*np.sinc(0*alpha) #Gleichanteil = a0
print("i = | ", end="")
for i in k:
    print("{0:>7d} | ".format(i), end="")
print("\nY = | ", end="")
for m in k:    
    print("{0:>7.2f} | ".format(A*alpha*np.sinc(m*alpha)), end="")
Y0 = A*alpha*np.sinc(k*alpha)

i = |       0 |       1 |       2 |       3 |       4 |       5 |       6 |       7 |       8 |       9 |      10 | 
Y = |  250.00 |  225.08 |  159.15 |   75.03 |    0.00 |  -45.02 |  -53.05 |  -32.15 |   -0.00 |   25.01 |   31.83 | 

2.2 Abschätzung des Spektrums durch FFT

Um das Spektrum des zeitkontinuierlichen Signals mit einer FFT abschätzen zu können, muss es zunächst abgetastet werden und zwar über eine ganzzahlige Anzahl Perioden $T_1$ (da es ja periodisch ist), wir starten zunächst mit einer Periode, $T_{mess1} = T_1$. Für eine effiziente Berechnung sollen $N_1 = 2^4 = 16$ Samples genommen werden.

VORBEREITUNG:

• Welche Abtastfrequenz $f_{S1}$ ist dafür erforderlich? $f_s = N1/T_1 = 16kHz$

• Welchen Abstand $T_{S1}$ haben die abgetasteten Punkte? In welchem Abstand liegen die Frequenzpunkte der FFT? $T_{S1} = \frac{1}{f_s} = 62.5 \, \mu\text{s}$, delta_f_fft = fs/N1 = 1kHz+ü

SIMULATION:

Zunächst einmal benötigen Sie einen passenden Zeitvektor z.B. mit t1 = np.arange(N_1)/f_S1 oder durch Abtasten des "zeitkontinuierlichen" Zeitvektors t (s.u.). Das abgetastete Signal $y_1[n]$ erzeugen Sie:

• Händisch mit y1 = np.ones(N_1); y1[4:N_1] = 0 oder durch y1 = np.concatenate((np.ones(4), np.zeros(12)))

• Durch Abtasten des "zeitkontinuierlichen" Signals ($f_S = 640 \text{ kHz}$) mit y1 = y[::40], mit diesem Befehl wird von y nur jedes vierzigste Element nach y1 kopiert. Das funktioniert natürlich nur, wenn wie hier $f_{S1} = f_S / 40$ ist. Wenn wie hier Start- und Endwert fehlen, werden alle Elemente von y vom ersten bis zum letzten berücksichtigt, Sie müssen ggf. noch Start- und Endwert anpassen, um den richtigen Ausschnitt finden.

• Durch Berechnung aus dem Zeitvektor mit Hilfe der np.where() Funktion wie im vorigen Lab.

Egal für welche Variante Sie sich entschieden haben, plotten Sie das abgetastete Signal mit dem folgenden Skript.

---

"""lab 1: fig, ax = plt.subplots(**size)
t = np.arange(0, 400e-6, T_S)
h = np.where(t<250e-6,1,0)
ax.stem(t,h, linefmt = "b-", markerfmt="ro", basefmt="k")
ax.set_xlabel(r"$t \; \rightarrow$ in s");"""

#Abtastung eines zeitkontinuierlichen Signals mithilfe einer FFT

fig, ax = plt.subplots(**size)
N1 = 16
f_S1 = 640e3/40
t1 = np.arange(N1)/f_S1 #step is T_S
y1 = np.concatenate((np.ones(4), np.zeros(12)))
#y1 = np.where(t1<20,1,0)
ax.stem(t1, y1, 'r')
ax.set_xlabel(r"$k T_S \; / \;\mathrm{s} \;\rightarrow$")
ax.set_ylabel(r"$y(k T_S)\; / \;\mathrm{mV} \;\rightarrow$");

<>:1: SyntaxWarning: invalid escape sequence '\;'
<>:1: SyntaxWarning: invalid escape sequence '\;'
C:\Users\taube\AppData\Local\Temp\ipykernel_30452\1562420728.py:1: SyntaxWarning: invalid escape sequence '\;'
  """lab 1: fig, ax = plt.subplots(**size)

2.2.1 FFT, erster Ansatz: $f_{S1} = 16 \text{ kHz}$, $N_1 = 16$

VORBEREITUNG:

• Welche höchste Frequenzkomponente ist im zeitkontinuierlichen Signal enthalten?

• Welche maximale Frequenzkomponente ist im abgetasteten Signal enthalten?

SIMULATION:

In numpy berechnen wir die FFT mit Y = np.fft.fft(y, N) (die FFT liefert komplexe Werte ...). Ohne den Parameter N wird die FFT über das gesamte Signal berechnet mit der Anzahl der Datenpunkte. Nutzen Sie das Notebook 03_DFT/DFT-Skalierung.ipynb zur korrekten Berechnung und Skalierung des Signals. Mit f = np.fft.fftfreq(N, T_S) erzeugen Sie eine passenden Frequenzvektor.

• Vervollständigen Sie das Skript in der nächsten Codezelle: Wie müssen Sie die Amplitude und die Frequenz skalieren, um physikalisch sinnvolle Werte zu erhalten?

fig, ax = plt.subplots(**size)
N1 = 16
f_S1 = N1 * 1e3
t1 = np.linspace(0.0, N1 / f_S1, num = N1, endpoint=False)
#y1 = ...
Y1 = np.fft.fft(y1, N1)/N1 * 1000
f1 = np.fft.fftfreq(N1, 1/f_S1)
ax.stem(f1, Y1, 'r')
ax.set_xlabel(r"$f \; / \;\mathrm{Hz} \;\rightarrow$")
ax.set_ylabel(r"$Y_1(F)\; / \;\mathrm{mV} \;\rightarrow$");
# ----------- Tabellenausgabe --------------
# Drucken Sie i, Y0 und Y1 aus

c:\Users\taube\anaconda3\envs\dsp_venv\Lib\site-packages\matplotlib\cbook.py:1699: ComplexWarning: Casting complex values to real discards the imaginary part
  return math.isfinite(val)
c:\Users\taube\anaconda3\envs\dsp_venv\Lib\site-packages\numpy\ma\core.py:3387: ComplexWarning: Casting complex values to real discards the imaginary part
  _data[indx] = dval
c:\Users\taube\anaconda3\envs\dsp_venv\Lib\site-packages\matplotlib\cbook.py:1345: ComplexWarning: Casting complex values to real discards the imaginary part
  return np.asarray(x, float)

Nach der SIMULATION:

• Vergleichen Sie die simulierten mit den berechneten Werten. Welche Koeffizienten weichen besonders stark von den erwarteten (berechneten) Werten ab? Was ist die Ursache dafür?

• Wie könnte man die Abweichungen reduzieren?

2.2.2 FFT mit Erhöhung der Abtastfrequenz auf $f_{S2} = 64 \text{ kHz}$ und $N_2$ Samples

Die Abtastfrequenz wird jetzt vervierfacht auf den Wert $f_{S2} = 64 \text{ kHz}$.

VORBEREITUNG:

• Wie lang muss das Messfenster $T_{Mess2}$ der FFT sein für eine Frequenzauflösung $\Delta f = 1$ kHz? Wieviele Samples $N_2$ sind dafür erforderlich? Antwort: $\Delta f = \frac{f_{S2}}{N_2} \text{daraus folgt:} \quad N_2 = \frac{f_{S2}}{\Delta f} = 64$ $T_{mess} = \frac {1}{\Delta f} = \frac {1}{1 \text{ kHz}} = 1 \text { ms}$

• Welche maximale Frequenzkomponente ist im abgetasteten Signal enthalten? Theoretisch fs/2 = 32 kHz, $f_{\text {max}} = \frac{N_2-1}{2*N2} = 31 \text { kHz} \text { weil reelwertig.. Symmetrie um} f_s/2$

SIMULATION:

• Kopieren Sie das Skript aus dem vorigen Unterpunkt und passsen Sie es an die ermittelten Parameter an.

• Vergleichen Sie erneut Rechnung mit Simulation

fig, ax = plt.subplots(**size)
N2 = 64
f_S2 = N2 * 1e3

t2 = 1/f_S2 
y2 = np.concatenate((np.ones(16), np.zeros(48)))
Y2 = np.fft.fft(y2, N2)/N2 *1000 # Umrechnung in mV
f2 = np.fft.fftfreq(N2, 1/f_S2)


ax.stem(f2, Y2, 'r')
ax.set_xlabel(r"$f \; / \;\mathrm{Hz} \;\rightarrow$")
ax.set_ylabel(r"$Y_1(F)\; / \;\mathrm{mV} \;\rightarrow$")

print("i    | ", end="")
for i in k:
    print("{0:>7d} | ".format(i), end="")

print("\nY0   | ", end="")
for i in Y0:
    print("{0:>7.2f} | ".format(abs(i)), end="")

print("\nY2   | ", end="")
for i in Y2:
    print("{0:>7.2f} | ".format(abs(i)), end="")

i    |       0 |       1 |       2 |       3 |       4 |       5 |       6 |       7 |       8 |       9 |      10 | 
Y0   |  250.00 |  225.08 |  159.15 |   75.03 |    0.00 |   45.02 |   53.05 |   32.15 |    0.00 |   25.01 |   31.83 | 
Y2   |  250.00 |  225.17 |  159.41 |   75.30 |    0.00 |   45.47 |   53.83 |   32.80 |    0.00 |   25.84 |   33.15 |   21.49 |    0.00 |   18.55 |   24.63 |   16.45 |    0.00 |   14.91 |   20.21 |   13.76 |    0.00 |   12.88 |   17.72 |   12.22 |    0.00 |   11.73 |   16.33 |   11.39 |    0.00 |   11.17 |   15.70 |   11.06 |    0.00 |   11.06 |   15.70 |   11.17 |    0.00 |   11.39 |   16.33 |   11.73 |    0.00 |   12.22 |   17.72 |   12.88 |    0.00 |   13.76 |   20.21 |   14.91 |    0.00 |   16.45 |   24.63 |   18.55 |    0.00 |   21.49 |   33.15 |   25.84 |    0.00 |   32.80 |   53.83 |   45.47 |    0.00 |   75.30 |  159.41 |  225.17 | 

2.2.3 FFT mit verlängertem Messfenster ($T_{Mess3} = 4 \text{ ms}$) und $f_{S3} = 16 \text{ kHz}$

Die Abtastfrequenz ist wieder wie zu Beginn $f_{S3} = 16 \text{ kHz}$, jetzt soll getestet werden wie sich eine Verlängerung des Messfensters um den Faktor 4 auswirkt.

VORBEREITUNG:

• Wie viele Samples $N_3$ benötigen Sie jetzt? Wieviele Perioden passen in Ihr Messfenster? $N_3 = f_s * T_{\text mess} = 16 \text{ kHz} * 4 \text{ ms} = 64$

• Welche Frequenzauflösung erzielt man mit diesem Setup? Welche Frequenzkomponente kann maximal dargestellt werden? $\Delta f = \frac{f_{S3}}{N_3} = \frac{16 \text{ kHz}}{64} = 250 \text{Hz}$ Theoretisch fs/2 = 8 kHz, $f_{\text {max}} = \frac{N_2-1}{2*N2} = 7,75 \text { kHz}$

SIMULATION:

Da Ihr Messfenster jetzt mehrere Perioden des Signals umfasst, müssen Sie Ihr Signal erzeugen entweder

• Händisch durch mehrfaches Aneinanderhängen eines Pulses z.B. mit y3 = np.tile(y1, 4)

• Durch Abtasten des "zeitkontinuierlichen" Signals wie am Anfang von 2.2 beschrieben. Achten Sie wieder darauf, Anfang und Ende passend zu wählen.

Plotten Sie Ihr Zeitsignal, damit Sie sicher sind dass Ihr Signal so aussieht wie Sie sich das vorstellen.

• Vergleichen Sie erneut Rechnung mit Simulation

• Ähnelt das Spektrum mehr dem ursprünglichen Spektrum $Y_1$ oder dem Spektrum mit erhöhter Abtastrate $Y_2$? Mehr Y2

• Welchen Vorteil hat die Verlängerung des Messfenster? Welchen Vorteil hat sie im hier vorliegenden Fall? Mehr Datenpunkte können erfasst werden, führt zu genauerer Abbildung

fig, ax = plt.subplots(**size)
N3 = 64
f_S3 = 16 * 1e3
t3 = np.arange(N3)/f_S3
y3 = np.tile(y1, 4)
ax.stem(t3, y3, 'r')
ax.set_xlabel(r"$k T_S \; / \;\mathrm{s} \;\rightarrow$")
ax.set_ylabel(r"$y(k T_S)\; / \;\mathrm{mV} \;\rightarrow$");

fig2, ax = plt.subplots(**size)
t3 = 1/f_S3 
Y3 = np.fft.fft(y3, N3)/N3 *1000 
f3 = np.fft.fftfreq(N3, 1/f_S3)


ax.stem(f3, Y3, 'r')
ax.set_xlabel(r"$f \; / \;\mathrm{Hz} \;\rightarrow$")
ax.set_ylabel(r"$Y_1(F)\; / \;\mathrm{mV} \;\rightarrow$")



# ----------- Tabellenausgabe --------------
# Drucken Sie i, Y0 und Y1 aus
# Tabellenausgabe

print("i    | ", end="")
for i in k:
    print("{0:>7d} | ".format(i), end="")

print("\nY0   | ", end="")
for i in Y0:
    print("{0:>7.2f} | ".format(abs(i)), end="")

print("\nY3   | ", end="")
for i in Y3:
    print("{0:>7.2f} | ".format(abs(i)), end="")

i    |       0 |       1 |       2 |       3 |       4 |       5 |       6 |       7 |       8 |       9 |      10 | 
Y0   |  250.00 |  225.08 |  159.15 |   75.03 |    0.00 |   45.02 |   53.05 |   32.15 |    0.00 |   25.01 |   31.83 | 
Y3   |  250.00 |    0.00 |    0.00 |    0.00 |  226.53 |    0.00 |    0.00 |    0.00 |  163.32 |    0.00 |    0.00 |    0.00 |   79.55 |    0.00 |    0.00 |    0.00 |    0.00 |    0.00 |    0.00 |    0.00 |   53.15 |    0.00 |    0.00 |    0.00 |   67.65 |    0.00 |    0.00 |    0.00 |   45.06 |    0.00 |    0.00 |    0.00 |    0.00 |    0.00 |    0.00 |    0.00 |   45.06 |    0.00 |    0.00 |    0.00 |   67.65 |    0.00 |    0.00 |    0.00 |   53.15 |    0.00 |    0.00 |    0.00 |    0.00 |    0.00 |    0.00 |    0.00 |   79.55 |    0.00 |    0.00 |    0.00 |  163.32 |    0.00 |    0.00 |    0.00 |  226.53 |    0.00 |    0.00 |    0.00 | 

2.2.4 FFT über nicht-ganzzahlige Anzahl von Perioden

Das Signal wird wieder mit $f_{S4} = 16\text{ kHz}$ abgetastet, allerdings ist das Messfenster jetzt nur noch $60 \, T_S$ lang.

VORBEREITUNG:

• Welche Länge $T_{Mess4}$ hat jetzt Ihr Messfenster und wieviele Signalperioden passen in das Messfenster? $T_{mess4} = 60 * T_S = 60 / f_S = 60 / 16\text{ kHz} = 3,75 ms$

• Welche Frequenzauflösung $\Delta f$ ergibt sich, in welchem Raster liegen die Frequenzbins $f_k$? $\Delta f = 266,67 \text { Hz}$

SIMULATION:

Signal und Zeitvektor erzeugen Sie am einfachsten durch Verkürzen der Signale aus dem letzten Abschnitt, also mit y4 = y3[:60] oder alternativ y4 = y3[:-4] (negative Indizes zählen vom Ende aus). Der Index, der das Ende markiert ("60" bzw. "-4") definiert dabei den ersten Wert, der nicht mit ausgegeben wird. Das ist z.B. bei arange(0,10,1) genauso und ist eine Folge des 0-based indexing (der erste Wert eines Arrays wird mit 0 angesprochen).

• Was hat sich im Vergleich zum Spektrum $Y_3$ des vorigen Abschnitts geändert? Leckeffekt, Lattenzauneffekt

fig, ax = plt.subplots(**size)
N4 = 60
f_S4 = 16 * 1e3
t4 = 1/f_S4
y4 = y3[:60]

Y4 = np.fft.fft(y4, N4)/N4 *1000 
f4 = np.fft.fftfreq(N4, 1/f_S4)


ax.stem(f4, Y4, 'r')
ax.set_xlabel(r"$f \; / \;\mathrm{Hz} \;\rightarrow$")
ax.set_ylabel(r"$Y_1(F)\; / \;\mathrm{mV} \;\rightarrow$")



# ----------- Tabellenausgabe --------------
# Drucken Sie i, Y0 und Y1 aus
# Tabellenausgabe

print("i    | ", end="")
for i in k:
    print("{0:>7d} | ".format(i), end="")

print("\nY0   | ", end="")
for i in Y0:
    print("{0:>7.2f} | ".format(abs(i)), end="")

print("\nY4   | ", end="")
for i in Y4:
    print("{0:>7.2f} | ".format(abs(i)), end="")

i    |       0 |       1 |       2 |       3 |       4 |       5 |       6 |       7 |       8 |       9 |      10 | 
Y0   |  250.00 |  225.08 |  159.15 |   75.03 |    0.00 |   45.02 |   53.05 |   32.15 |    0.00 |   25.01 |   31.83 | 
Y4   |  266.67 |   18.52 |   26.52 |   62.62 |  212.93 |   55.77 |   51.29 |  113.09 |   99.64 |   34.91 |   28.87 |   81.28 |   16.67 |    4.41 |    1.45 |    0.00 |    1.30 |    3.57 |   12.11 |   52.79 |   16.67 |   17.79 |   44.36 |   43.41 |   16.67 |   14.94 |   45.26 |    9.92 |    2.79 |    0.97 |    0.00 |    0.97 |    2.79 |    9.92 |   45.26 |   14.94 |   16.67 |   43.41 |   44.36 |   17.79 |   16.67 |   52.79 |   12.11 |    3.57 |    1.30 |    0.00 |    1.45 |    4.41 |   16.67 |   81.28 |   28.87 |   34.91 |   99.64 |  113.09 |   51.29 |   55.77 |  212.93 |   62.62 |   26.52 |   18.52 | 

3. Schätzung des Fourierspektrums von Sinustönen mit überlagertem Rauschen durch FFT

In diesem Versuchsteil versuchen wir mit Hilfe der FFT Amplituden und Frequenzen von zwei Sinustönen zu schätzen mit $f_1 = 990\text{ Hz}$ und $A_1 = 100\text{ mV}$ sowie $f_2 = 1010\text{ Hz}$ und $A_2 = 2\text{ mV}$, überlagert ist eine gleichverteilte Rauschstörung im Bereich $\pm 1 \text{ mV}$.

Tipp: Diesen Versuchsteil können Sie auch mit pyfda durchführen, im Tab y[n] wählen Sie dazu ein Sinussignal mit überlagertem Rauschen aus und schauen sich das Spektrum im Untertab "Frequency" an. Enablen Sie "Stimulus X" und deaktivieren Sie "Response Y".

3.1 Abtastfrequenz frei wählbar

Die Anzahl der Datenpunkte $N$ für die FFT und die Abtastfrequenz $f_S$ können zunächst frei gewählt werden.

VORBEREITUNG:

• Welche Abtastfrequenz benötigen Sie mindestens?

  Die Abtastfrequenz muss größer als $2 f_2 = 2020 \text{ Hz}$ sein.

• Wählen Sie die kleinstmöglichen Werte für $N$ und $f_S$ so, dass kein Frequenzfehler und kein Leckeffekt auftritt. Tipp: Bestimmen Sie mit $T_{min} = L_1 T_1 = L_2 T_2$ zunächst das kürzeste Zeitfenster, in dem sowohl $f_1$ als auch $f_2$ periodisch sind. Danach können Sie mit $T_{Mess} = \frac{N}{f_S} = M T_{min}$ eine Bedingung für $N$ und $T_S$ ableiten $$T_{\text{min}} = L_1 * T_1 = L_2 * T_2 \text{ daraus folgt: } \frac{L_1}{L_2} = \frac{T_2}{T_1} = \frac{99}{101} $$ $$T_{\text{min}} = 99 \times T_1 = 100 \text{ ms} \text{ und } f_{\text{s,min}} = \frac{f_s}{2} + \Delta f = \frac{f_s}{2} * \frac{1}{T_{\text{mess}}} = 2030 \text{ Hz} \text{ folgt: } $$ $$\frac{N}{f_{\text{s,min}}} \stackrel{!}{=} 100 \text{ ms, folglich: } N_{\text{min}} = 203 $$

• Wie groß ist das Verhältnis der Amplituden und wie groß das Verhältnis der Leistungen beider Sinustöne in dB?

SIMULATION:

• Wählen Sie $f_S = 2500 \text{ Hz}$ und $N = 250$ für ein leichter ablesbares Spektrum.

• Stellen Sie das Spektrum als (Amplituden-)Betragsspektrum im logarithmischen Maßstab dar, die Werte sollen relativ zum Maximum dargestellt werden. Zur Verbesserung der Darstellung können Sie mit np.fft.fftshift positive und negative Achse von Frequenzvektor und Spektrum vertauschen.

N5 =  250
f_S5 = 2500
t = np.arange(0,0.1,0.1/250)
noi = 2 * np.random.rand(N5) - 1
y5 = 100 * np.sin(2 * np.pi * 990 * t) + 2 * np.sin(2 * np.pi * 1010 * t) + noi
t5 = 1/f_S5


fig, ax = plt.subplots(**size)
Y5 = 20*np.log((np.fft.fft(y5/100,N5)/N5))
f5 = np.fft.fftfreq(N5, t5)
ax.stem(f5, Y5, 'r')
ax.set_xlabel(r"$f \; / \;\mathrm{Hz} \;\rightarrow$")
ax.set_ylabel(r"$Y_5(F)\; / \;\mathrm{dB} \;\rightarrow$");

3.2 Feste Abtastfrequenz $f_S = 4 \text{ kHz}$

Die Abtastfrequenz ist jetzt fest eingestellt auf $f_S = 4 \text{ kHz}$, die Anzahl der Datenpunkte $N$ für die FFT soll eine Zweierpotenz sein, $N = 2^L$, für eine möglichst effiziente Berechnung der FFT.

Auch diesen Unterpunkt können Sie mit pyfda bearbeiten.

VORBEREITUNG:

• Gibt es bei der gewählten Abtastfrequenz einen Wert für $N$, bei dem kein Frequenzfehler und kein Leckeffekt auftritt (mit und ohne Beschränkung auf Zweierpotenzen)? Gibt es für $N = 256$ eine Abtastfrequenz, bei der kein Leckeffekt auftritt? $$N = T_{\text mess} * f_S = 400$$

$$f_S = \frac{N}{T_{\text mess}} = 2,56 \text{ kHz}$$

SIMULATION:

• Der Leckeffekt soll jetzt durch Verwendung "weicherer" Fenster reduziert werden. Fensterfunktionen finden Sie unter scipy.signal.windows (https://docs.scipy.org/doc/scipy/reference/signal.windows.html). Vor Berechnung der FFT müssen Sie dazu Ihr Signal mit der Fensterfunktion multiplizieren, also z.B. y_win = y * sig.windows.hann(256, sym=False). Dabei wurde angenommen, dass auch Ihr Zeitsignal y eine Länge von $N = 256$ hat. Für die Spektralanalyse periodischer Signale (im Gegensatz zum Filterdesign und zur Analyse von impulsförmigen Signalen) muss sym = False gewählt werden.

• Testen Sie den Einfluss der Fensterlänge $N = 256$, $N=512$ und $N=1024$ ... und den Effekt eines rechteckförmigen (boxcar bzw. gar keine Fensterung) und eines Hann-Fensters. Wie genau werden die Amplituden und Frequenzen geschätzt? Schauen Sie sich auch die anderen Fensterfunktionen an.

N6 =  256
f_S6 = 4000
t = np.arange(0, N6/f_S6, 1/f_S6)
noi = 2 * np.random.rand(N6) - 1
y6 = 100 * np.sin(2 * np.pi * 990 * t) + 2 * np.sin(2 * np.pi * 1010 * t) + noi
y_win = y6 * sig.windows.hann(N6, sym=False)
y_boxcar = y6*sig.windows.boxcar(N6, sym = False)

print("Mit N=256")

fig, ax = plt.subplots(**size)
t6 = 1/f_S6
Y6 = (np.fft.fft(y6,N6)/N6)
f6 = np.fft.fftfreq(N6, t6)
ax.stem(f6, Y6, 'r')
ax.set_xlabel(r"$f \; / \;\mathrm{Hz} \;\rightarrow$")
ax.set_ylabel(r"$Y_6(F)\; / \;\mathrm{mV} \;\rightarrow$");
ax.set_title("Ohne Filter N=256")


fig, ax = plt.subplots(**size)
Y6 = (np.fft.fft(y_win,N6)/N6)
f6 = np.fft.fftfreq(N6, t6)
ax.stem(f6, Y6, 'r')
ax.set_xlabel(r"$f \; / \;\mathrm{Hz} \;\rightarrow$")
ax.set_ylabel(r"$Y_6(F)\; / \;\mathrm{mV} \;\rightarrow$");
ax.set_title("Hann N=256")

fig, ax = plt.subplots(**size)
Y6 = (np.fft.fft(y_boxcar,N6)/N6)
f6 = np.fft.fftfreq(N6, t6)
ax.stem(f6, Y6, 'r')
ax.set_xlabel(r"$f \; / \;\mathrm{Hz} \;\rightarrow$")
ax.set_ylabel(r"$Y_6(F)\; / \;\mathrm{mV} \;\rightarrow$");
ax.set_title("Boxcar N=256")


N6 =  512
f_S6 = 4000
t = np.arange(0, N6/f_S6, 1/f_S6)
noi = 2 * np.random.rand(N6) - 1
y6 = 100 * np.sin(2 * np.pi * 990 * t) + 2 * np.sin(2 * np.pi * 1010 * t) + noi
y_win = y6 * sig.windows.hann(N6, sym=False)
y_boxcar = y6*sig.windows.boxcar(N6, sym = False)

print("Mit N=512")

fig, ax = plt.subplots(**size)
t6 = 1/f_S6
Y6 = (np.fft.fft(y6,N6)/N6)
f6 = np.fft.fftfreq(N6, t6)
ax.stem(f6, Y6, 'r')
ax.set_xlabel(r"$f \; / \;\mathrm{Hz} \;\rightarrow$")
ax.set_ylabel(r"$Y_6(F)\; / \;\mathrm{mV} \;\rightarrow$");
ax.set_title("Ohne Filter N=512")


fig, ax = plt.subplots(**size)
Y6 = (np.fft.fft(y_win,N6)/N6)
f6 = np.fft.fftfreq(N6, t6)
ax.stem(f6, Y6, 'r')
ax.set_xlabel(r"$f \; / \;\mathrm{Hz} \;\rightarrow$")
ax.set_ylabel(r"$Y_6(F)\; / \;\mathrm{mV} \;\rightarrow$");
ax.set_title("Hann N=512")

fig, ax = plt.subplots(**size)
Y6 = (np.fft.fft(y_boxcar,N6)/N6)
f6 = np.fft.fftfreq(N6, t6)
ax.stem(f6, Y6, 'r')
ax.set_xlabel(r"$f \; / \;\mathrm{Hz} \;\rightarrow$")
ax.set_ylabel(r"$Y_6(F)\; / \;\mathrm{mV} \;\rightarrow$");
ax.set_title("Boxcar N=512")

N6 =  1024
f_S6 = 4000
t = np.arange(0, N6/f_S6, 1/f_S6)
noi = 2 * np.random.rand(N6) - 1
y6 = 100 * np.sin(2 * np.pi * 990 * t) + 2 * np.sin(2 * np.pi * 1010 * t) + noi
y_win = y6 * sig.windows.hann(N6, sym=False)
y_boxcar = y6*sig.windows.boxcar(N6, sym = False)

print("Mit N=1024")

fig, ax = plt.subplots(**size)
t6 = 1/f_S6
Y6 = (np.fft.fft(y6,N6)/N6)
f6 = np.fft.fftfreq(N6, t6)
ax.stem(f6, Y6, 'r')
ax.set_xlabel(r"$f \; / \;\mathrm{Hz} \;\rightarrow$")
ax.set_ylabel(r"$Y_6(F)\; / \;\mathrm{mV} \;\rightarrow$");
ax.set_title("Ohne Filter N=1024")


fig, ax = plt.subplots(**size)
Y6 = (np.fft.fft(y_win,N6)/N6)
f6 = np.fft.fftfreq(N6, t6)
ax.stem(f6, Y6, 'r')
ax.set_xlabel(r"$f \; / \;\mathrm{Hz} \;\rightarrow$")
ax.set_ylabel(r"$Y_6(F)\; / \;\mathrm{mV} \;\rightarrow$");
ax.set_title("Hann N=1024")

fig, ax = plt.subplots(**size)
Y6 = (np.fft.fft(y_boxcar,N6)/N6)
f6 = np.fft.fftfreq(N6, t6)
ax.stem(f6, Y6, 'r')
ax.set_xlabel(r"$f \; / \;\mathrm{Hz} \;\rightarrow$")
ax.set_ylabel(r"$Y_6(F)\; / \;\mathrm{mV} \;\rightarrow$");
ax.set_title("Boxcar N=1024")

Mit N=256
Mit N=512
Mit N=1024

Text(0.5, 1.0, 'Boxcar N=1024')

4. Schätzung des Fourierspektrums eines rechteckförmigen Impulses durch FFT

In diesem Versuchsteil schätzen wir das Fourierintegral des rechteckförmigen Impulses $y(t)$ mit $\Delta T = 250 \,\mu$s und $A = 1000\,\text{mV}$ (s. nächster Plot) mit Hilfe der FFT.

4.1 Berechnung der Fouriertransformierten

VORBEREITUNG:

• Berechnen Sie die Fouriertransformierte des Impulses und skizzieren Sie die Betragsfunktion (oder plotten Sie sie). $$U(f) = 250 \mu \text{V} \times \text{si}(\pi \times 250 \mu s \times f)$$

• Bei welcher Frequenz hat die Fouriertransformierte die erste Nullstelle?

  Die erste Nullstelle $f_0$ der Fouriertransformierten ist bei $\pi f_0 \Delta T_1 = \pi \;\Rightarrow \; f_0 = 1/ \Delta T = 4\, \text{kHz}$.

fig, ax = plt.subplots(**size)
Delta_T = 250e-6
A = 1000
t = np.arange(-200e-6, 200e-6, 1e-6)
ax.plot(t*1e6, A*np.where((t >= -Delta_T/2) & (t < Delta_T/2), 1, 0), 'r', lw=2)
ax.axvline(0, ls='-',color='k')
ax.set_xlabel(r"$t \; / \;\mu\mathrm{s} \;\rightarrow$")
ax.set_ylabel(r"$y(t) \; \mathrm{in \;mV} \;\rightarrow$");

Das zeitkontinuierliche Signal $y(t)$ hat eine endliche zeitliche Ausdehnung und endliche Energie, daher kann sein Amplitudenspektrum mit dem Fourierintegral beschrieben werden:

$$Y\left( f \right) = \int_{-\infty}^{\infty}y(t)\text{e}^{-\text{j} 2 \pi f t} \text {d}t$$

Der Impuls ist symmetrisch zur y-Achse was die Berechnung deutlich vereinfacht:

$$ \begin{align} Y\left( f \right) &= A \int_{-\Delta T/2}^{\Delta T/2}\text{e}^{-\text{j} 2 \pi f t} \text {d}t = \left. \frac A{-\text{j} 2 \pi f} \text{e}^{-\text{j} 2 \pi f t}\right|_{-\Delta T/2} ^{\Delta T/2} = A\frac{\text{e}^{\text{j} \pi f \Delta T}- \text{e}^{-\text{j} \pi f \Delta T}}{2 \pi \text{j}f} = A\frac{\sin \pi f \Delta T}{\pi f} \\ &= A\Delta T \frac{\sin \pi f \Delta T}{\pi f \Delta T} = A\Delta T \text{sinc} f \Delta T \text{ mit } \text{sinc} f = \frac{\sin \pi f}{\pi f} \end{align}$$ Das Ergebnis ist die Amplitudenspektrumsdichte in V/Hz.

Berechnet man die FFT anstatt des Fourierintegrals, überstreicht jeder Frequenzpunkt die Bandbreite $\Delta f = f_S/N$. Möchte man die FFT eines Energiesignals so skalieren dass man (in etwa) identische Zahlenwerte erhält wie bei der spektralen Amplitudendichte des ursprünglichen zeitkontinuierlichen Signals, so muss man mit $1/\Delta f = N/f_S$ multiplizieren. Da die FFT bereits mit $1/N$ skaliert ist, hebt sich der Faktor $N$ auf. Man muss also nur durch $f_S$ dividieren.

Für unendlich ausgedehnte Leistungssignale (periodische Signale, stationäre und nicht-stationäre Prozesse) konvergiert das Integral nicht, man nimmt stattdessen die Fourierreihe (periodische Signale) oder die spektrale Leistungsdichte, die über die Autokorrelationsfunktion berechnet wird (aber nicht hier und jetzt ...).

fig, ax = plt.subplots(**size)
Delta_T = 250e-6
A = 1000
f = np.arange(-2.2e4, 2.2e4, 1e2)
ax.plot(f, A * Delta_T*np.sinc(f * Delta_T), 'r', lw=2)
ax.axvline(0, ls='-',color='k'); ax.axhline(0, ls='-',color='k') 
ax.set_title(r"Fourierintegral des Rechteckimpulses mit $\Delta T = 250\, \mu \mathrm{s}$")
ax.set_xlabel(r"$f \; / \;\mathrm{Hz} \;\rightarrow$")
ax.set_ylabel(r"$Y(f) \; \mathrm{in \;mV/Hz} \;\rightarrow$");

4.2 Schätzung der Fouriertransformierten mit Hilfe der FFT

Im folgenden soll die Fouriertransformierte des Impulses mit Hilfe einer FFT abgeschätzt werden, dazu wird der Impuls mit $L = 16$ Samples abgetastet.

VORBEREITUNG:

• Welche Abtastfrequenz $f_S$ ist dazu notwendig?

  $f_S = 16/T_{\text mess} = 40 \text{ kHz}$

• Wie groß ist die Frequenzauflösung $\Delta f$? Wieviele Frequenzpunkte erhält man zwischen $f=0$ und der ersten Nullstelle $f_0$?

  $\Delta f = f_s/16 = 2500 Hz$

• Kann man die graphische Darstellung des Amplitudenspektrums zwischen $f= 0$ und der ersten Nullstelle $f_0$ durch Änderung der Abtastfrequenz verbessern? Nein

• Wie kann die graphische Darstellung des Betragsgangs verbessert werden ohne die Abtastfrequenz zu verändern? Erhöhung der Anzahl der Datenpunkte

• Wie müssen die Ergebnisse skaliert werden, um physikalisch korrekte Werte für die spektrale Amplitudendichte in V/Hz zu bekommen (so wie theoretisch berechnet)? Skalieren mit 1/f_s

  SIMULATION:

• Berechnen Sie die FFT des Rechteckpulses mit $L=16$ und mit $N_{FFT} = 2^9 = 512$ Punkten.

• Erstellen Sie einen Plot mit der FFT mit $L=16$ Punkten, $N_{FFT} = 512$ Punkten und dem idealen Amplitudenspektrum des Rechteckimpulses.

fig, ax = plt.subplots(**size)
f_S = 64e3
Delta_T = 250e-6
A = 1000
f = np.arange(-32e3, 32e3, 1e2)
f_16 = np.fft.fftshift(np.fft.fftfreq(16, d= 1/f_S))
f_512 = np.fft.fftshift(np.fft.fftfreq(512, d= 1/f_S))
y_id = A * Delta_T*np.sinc(f * Delta_T)
y_16 = A*np.ones(16)
Y_16 = np.fft.fftshift(np.abs(np.fft.fft(y_16))/f_S)
Y_512 = np.fft.fftshift(np.abs(np.fft.fft(y_16, 512))/f_S)
ax.plot(f, y_id, 'r', lw=2, label="$Y_{id}$")
ax.plot(f_16, Y_16, 'bo', label="$Y_{16}$")
ax.plot(f_512, Y_512, label="$Y_{512}$")
ax.legend()
ax.axvline(0, ls='-',color='k'); ax.axhline(0, ls='-',color='k') 
ax.set_xlabel(r"$f \; / \;\mathrm{Hz} \;\rightarrow$")
ax.set_ylabel(r"$Y(f) \; \mathrm{in \;V/Hz} \;\rightarrow$");

---

Copyright

(c) 2016 - 2021 Prof. Dr. Christian Münker

This jupyter notebook is part of a collection of notebooks on various topics of Digital Signal Processing. The latest version can be found at https://github.com/chipmuenk/dsp.

This notebook is provided as Open Educational Resource. Feel free to use it for your own purposes. The text is licensed under Creative Commons Attribution 4.0, the code of the IPython examples under the MIT license. Please attribute the work as follows: Christian Münker, Digital Signal Processing - Vorlesungsunterlagen mit Simulationsbeispielen, 2020.